Load the necessary libraries
library(rstanarm) #for fitting models in STAN
library(brms) #for fitting models in STAN
library(coda) #for diagnostics
library(bayesplot) #for diagnostics
library(rstan) #for interfacing with STAN
library(emmeans) #for marginal means etc
library(broom) #for tidying outputs
library(DHARMa) #for residual diagnostics
library(tidybayes) #for more tidying outputs
library(ggeffects) #for partial plots
library(tidyverse) #for data wrangling etc
To investigate differential metabolic plasticity in barramundi (Lates calcarifer), Norin, Malte, and Clark (2015) exposed juvenile barramundi to various environmental changes (increased temperature, decreased salinity and increased hypoxia) as well as control conditions. Metabolic plasticity was calculated as the percentage difference in standard metabolic rate between the various treatment conditions and the standard metabolic rate under control conditions. They were interested in whether there was a relationship between metabolic plasticity and typical (control) metabolism and how the different treatment conditions impact on this relationship.
A total of 60 barramundi juveniles were subject to each of the three conditions (high temperature, low salinity and hypoxia) in addition to control conditions. Fish mass was also recorded as a covariate as this is known to influence metabolic parameters.
Barramundi
Format of norin.csv data files
| FISHID | MASS | TRIAL | SMR_contr | CHANGE |
|---|---|---|---|---|
| 1 | 35.69 | LowSalinity | 5.85 | -31.92 |
| 2 | 33.84 | LowSalinity | 6.53 | 2.52 |
| 3 | 37.78 | LowSalinity | 5.66 | -6.28 |
| .. | .. | .. | .. | .. |
| 1 | 36.80 | HighTemperature | 5.85 | 18.32 |
| 2 | 34.98 | HighTemperature | 6.53 | 19.06 |
| 3 | 38.38 | HighTemperature | 5.66 | 19.03 |
| .. | .. | .. | .. | .. |
| 1 | 45.06 | Hypoxia | 5.85 | -18.61 |
| 2 | 43.51 | Hypoxia | 6.53 | -5.37 |
| 3 | 45.11 | Hypoxia | 5.66 | -13.95 |
| FISHID | Categorical listing of the individual fish that are repeatedly sampled |
| MASS | Mass (g) of barramundi. Covariate in analysis |
| TRIAL | Categorical listing of the trial (LowSalinity: 10ppt salinity; HighTemperature: 35 degrees; Hypoxia: 45% air-sat. oxygen. |
| SMR_contr | Standard metabolic rate (mg/h/39.4 g of fish) under control trial conditions (35 ppt salinity, 29 degrees, normoxia) |
| CHANGE | Percentage difference in Standard metabolic rate (mg/h/39.4 g of fish) between Trial conditions and control adjusted for 'regression to the mean'. |
norin = read_csv('../public/data/norin.csv', trim_ws=TRUE)
##
## ── Column specification ───────────────────────
## cols(
## FISHID = col_double(),
## MASS = col_double(),
## TRIAL = col_character(),
## SMR_contr = col_double(),
## CHANGE = col_double()
## )
glimpse(norin)
## Rows: 180
## Columns: 5
## $ FISHID <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1…
## $ MASS <dbl> 35.69, 33.84, 37.78, 26.58, 37.62, 37.68, 30.62, 50.37, 24.9…
## $ TRIAL <chr> "LowSalinity", "LowSalinity", "LowSalinity", "LowSalinity", …
## $ SMR_contr <dbl> 5.847466, 6.530707, 5.659556, 6.278200, 4.407336, 4.818589, …
## $ CHANGE <dbl> -31.919389, 2.520929, -6.284968, -4.346675, -3.071329, -15.0…
Model formula: \[ y_i \sim{} \mathcal{N}(\mu_i, \sigma^2)\\ \mu_i =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]
where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of temperature and (centered) mean fish size on SDA peak. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual fish.
norin = norin %>% mutate(FISHID=factor(FISHID),
TRIAL=factor(TRIAL))
ggplot(norin, aes(y=CHANGE, x=TRIAL)) + geom_boxplot()
ggplot(norin, aes(y=CHANGE, x=SMR_contr, shape=TRIAL, color=TRIAL)) +
geom_smooth(method='lm') + geom_point()
#ggplot(norin, aes(y=CHANGE, x=MASS, shape=TRIAL, color=TRIAL)) +
#geom_smooth(method='lm') + geom_point()
ggplot(norin, aes(y=CHANGE, x=as.numeric(FISHID), color=TRIAL)) +
geom_point() + geom_line()
#ggplot(norin, aes(y=MASS, x=TRIAL)) + geom_boxplot()
ggplot(norin, aes(y=CHANGE, x=MASS, color=TRIAL)) + geom_point() + geom_smooth(method='lm')
norin.rstanarm = stan_glmer(CHANGE ~ (1|FISHID)+TRIAL*SMR_contr+MASS, data=norin,
prior_PD=TRUE,
iter=5000, warmup=2000, chains=3, thin=5, refresh=0)
prior_summary(norin.rstanarm)
posterior_vs_prior(norin.rstanarm, color_by='vs', group_by=TRUE,
facet_args=list(scales='free_y'), pars=c('(Intercept)'))
ggpredict(norin.rstanarm, ~TRIAL*SMR_contr) %>% plot(add.data=TRUE)
norin.rstanarm %>% get_variables()
plot(norin.rstanarm, 'mcmc_trace', regex_pars='^.Intercept|TRIAL|SMR|MASS|[sS]igma')
plot(norin.rstanarm, 'mcmc_acf_bar', regex_pars='^.Intercept|TRIAL|SMR|MASS|[sS]igma')
plot(norin.rstanarm, 'mcmc_rhat_hist', regex_pars='^.Intercept|TRIAL|SMR|MASS|[sS]igma')
plot(norin.rstanarm, 'mcmc_neff_hist', regex_pars='^.Intercept|TRIAL|SMR|MASS|[sS]igma')
#norin.rstan1 = stan_glmer(CHANGE ~ (TRIAL|FISHID)+TRIAL*SMR_contr+MASS, data=norin,
# iter=5000, warmup=2000, chains=3, thin=5, refresh=0, cores=3)
norin.rstanarm1 = stan_glmer(CHANGE ~ (SMR_contr|FISHID) + TRIAL*SMR_contr+MASS, data=norin,
prior_PD=FALSE,
iter=5000, warmup=2000, chains=3, thin=5, refresh=0, cores=3)
norin.rstanarm1 = update(norin.rstanarm1, prior_PD=FALSE)
norin.rstanarm2 = stan_glmer(CHANGE ~ (TRIAL*SMR_contr|FISHID) + TRIAL*SMR_contr+MASS, data=norin,
prior_PD=FALSE,
iter=5000, warmup=2000, chains=3, thin=5, refresh=0, cores=3)
posterior_vs_prior(norin.rstanarm1, color_by='vs', group_by=TRUE,
facet_args=list(scales='free_y'), pars=c('(Intercept)'))
ggpredict(norin.rstanarm1, ~TRIAL*SMR_contr) %>% plot(add.data=TRUE)
norin.rstanarm1 %>% get_variables()
plot(norin.rstanarm1, 'mcmc_trace', regex_pars='^.Intercept|TRIAL|^SMR|MASS|[sS]igma')
plot(norin.rstanarm1, 'mcmc_acf_bar', regex_pars='^.Intercept|TRIAL|^SMR|MASS|[sS]igma')
plot(norin.rstanarm1, 'mcmc_rhat_hist', regex_pars='^.Intercept|TRIAL|^SMR|MASS|[sS]igma')
plot(norin.rstanarm1, 'mcmc_neff_hist', regex_pars='^.Intercept|TRIAL|^SMR|MASS|[sS]igma')
(l.1 <- loo(norin.rstanarm))
(l.2 <- loo(norin.rstanarm1))
loo_compare(l.1, l.2)
preds <- posterior_predict(norin.rstanarm, nsamples=250, summary=FALSE)
norin.resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = norin$CHANGE,
fittedPredictedResponse = apply(preds, 2, median))
plot(norin.resids)
g=ggpredict(norin.rstanarm) %>% plot
do.call('grid.arrange', g)
#ggemmeans(norin.rstan, ~TRIAL)
summary(norin.rstanarm)
nms <- norin.rstanarm1 %>% get_variables()
wch <- grep('^.Intercept|TRIAL|^SMR|[sS]igma', nms)
tidyMCMC(norin.rstanarm$stanfit,conf.int=TRUE, conf.method='HPDinterval',
rhat=TRUE, ess=TRUE, pars=nms[wch], estimate.method='median')
tidyMCMC(norin.rstanarm1$stanfit,conf.int=TRUE, conf.method='HPDinterval',
rhat=TRUE, ess=TRUE, pars=nms[wch], estimate.method='median')
norin.grid = with(norin, list(SMR_contr=seq(min(SMR_contr),max(SMR_contr), len=100)))
newdata = emmeans(norin.rstanarm, ~TRIAL|SMR_contr, at=norin.grid) %>% as.data.frame
head(newdata)
ggplot(newdata, aes(y=emmean, x=SMR_contr, color=TRIAL)) +
geom_ribbon(aes(ymin=lower.HPD, ymax=upper.HPD, fill=TRIAL), alpha=0.3,color=NA) +
geom_line()
norin.grid = with(norin, list(SMR_contr=c(min(SMR_contr),mean(SMR_contr),max(SMR_contr))))
emmeans(norin.rstan, pairwise~TRIAL|SMR_contr, at=norin.grid)
norin.em = emmeans(norin.rstan, pairwise~TRIAL|SMR_contr, at=norin.grid)$contrast %>%
gather_emmeans_draws() %>%
mutate(Fit=.value)
norin.em
norin.em %>% group_by(contrast) %>% median_hdci(Fit)
norin.em %>% group_by(contrast, SMR_contr) %>% median_hdci(Fit)
## norin.em %>%
## group_by(contrast) %>%
## summarize(P=sum(Fit>0)/n())
norin.em %>%
group_by(contrast, SMR_contr) %>%
summarize(P=sum(Fit>0)/n())
bayes_R2(norin.rstanarm, re.form=NA) %>% median_hdi
bayes_R2(norin.rstanarm, re.form=~(1|FISHID)) %>% median_hdi
#bayes_R2(norin.rstan1, re.form=~(SMR_contr|FISHID)) %>% median_hdi
norin = norin %>% mutate(FISHID=factor(FISHID),
TRIAL=factor(TRIAL))
ggplot(norin, aes(y=CHANGE, x=TRIAL)) + geom_boxplot()
ggplot(norin, aes(y=CHANGE, x=SMR_contr, shape=TRIAL, color=TRIAL)) +
geom_smooth(method='lm') + geom_point()
ggplot(norin, aes(y=CHANGE, x=MASS, shape=TRIAL, color=TRIAL)) +
geom_smooth(method='lm') + geom_point()
ggplot(norin, aes(y=CHANGE, x=as.numeric(FISHID), color=TRIAL)) +
geom_point() + geom_line()
##ggplot(norin, aes(y=MASS, x=TRIAL)) + geom_boxplot()
##ggplot(norin, aes(y=CHANGE, x=MASS, color=TRIAL)) + geom_point() + geom_smooth(method='lm')
norin %>% group_by(TRIAL) %>%
summarise(median(CHANGE),
mad(CHANGE))
priors <- prior(normal(50, 20), class='Intercept') +
prior(normal(0, 10), class='b') +
prior(gamma(2,1), class='sigma') +
prior(gamma(2,1), class='sd')
norin.form <- bf(CHANGE ~ (1|FISHID)+TRIAL*SMR_contr+MASS,
family=gaussian)
norin.brm1 = brm(norin.form,
data=norin,
prior = priors,
sample_prior = 'yes',
iter=5000, warmup=2000,
chains=3, thin=5, refresh=0)
norin.brm1 %>% get_variables()
pars <- norin.brm1 %>% get_variables()
wch <- grepl('^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd', pars, perl=TRUE)
stan_trace(norin.brm1$fit, pars=pars[wch])
stan_ac(norin.brm1$fit, pars=pars[wch])
stan_rhat(norin.brm1$fit, pars=pars[wch])
stan_ess(norin.brm1$fit, pars=pars[wch])
##mcmc_plot(norin.brms, type='trace',
## regex_pars='^b.*|sigma|^sd')
##mcmc_plot(norin.brms, type='trace',
## regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
##mcmc_plot(norin.brms, type='acf_bar',
## regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
##mcmc_plot(norin.brms, type='rhat_hist',
## regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
##mcmc_plot(norin.brms, type='neff_hist',
## regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
preds <- posterior_predict(norin.brm1, nsamples=250, summary=FALSE)
norin.resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = norin$CHANGE,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse =FALSE)
plot(norin.resids)
#norin.rstan1 = stan_glmer(CHANGE ~ (TRIAL|FISHID)+TRIAL*SMR_contr+MASS, data=norin,
# iter=5000, warmup=2000, chains=3, thin=5, refresh=0, cores=3)
norin.form <- bf(CHANGE ~ (TRIAL|FISHID) + TRIAL*SMR_contr+MASS,
family=gaussian)
norin.brm2 = brm(norin.form, data=norin,
prior = priors,
sample_prior = 'yes',
iter=5000, warmup=2000,
chains=3, thin=5, refresh=0, cores=3,
control=list(adapt_delta=0.99))
norin.brm2 %>% get_variables()
pars <- norin.brm2 %>% get_variables()
## wch <- grepl('^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd', pars, perl=TRUE)
wch <- grepl('^b_.*|[sS]igma|^sd_.*', pars, perl=TRUE)
stan_trace(norin.brm2$fit, pars=pars[wch])
stan_ac(norin.brm2$fit, pars=pars[wch])
stan_rhat(norin.brm2$fit)#, pars=pars[wch])
stan_ess(norin.brm2$fit)#, pars=pars[wch])
##mcmc_plot(norin.brm2, type='trace',
## regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
##mcmc_plot(norin.brm2, type='trace',
## regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
##mcmc_plot(norin.brm2, type='acf_bar',
## regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
##mcmc_plot(norin.brm2, type='rhat_hist',
## regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
##mcmc_plot(norin.brm2, type='neff_hist',
## regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
(l.1 <- loo(norin.brm1))
(l.2 <- loo(norin.brm2))
loo_compare(l.1, l.2)
preds <- posterior_predict(norin.brm2, nsamples=250, summary=FALSE)
norin.resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = norin$CHANGE,
fittedPredictedResponse = apply(preds, 2, median))
plot(norin.resids)
g <- norin.brm2 %>%
conditional_effects() %>%
plot(points=TRUE, ask=FALSE)
library(patchwork)
g[[1]] + g[[2]] + g[[3]] + g[[4]]
##g=ggpredict(norin.brms1) %>% plot
##library(patchwork)
##g[[1]] + g[[2]] + g[[3]]
##do.call('grid.arrange', g)
ggemmeans(norin.brm2, ~TRIAL) %>% plot
summary(norin.brm2)
tidyMCMC(norin.brm2$fit,conf.int=TRUE, conf.method='HPDinterval',
rhat=TRUE, ess=TRUE, estimate.method='median') %>%
slice(1:11)
pars <- norin.brm2 %>% get_variables()
wch <- grep('^b.Intercept|TRIAL|^b.*SMR|[sS]igma|^sd', pars)
tidyMCMC(norin.brms1$fit,conf.int=TRUE, conf.method='HPDinterval',
rhat=TRUE, ess=TRUE, pars=pars[wch], estimate.method='median')
bayes_R2(norin.brm2, re.form=NA, summary=FALSE) %>%
median_hdci
bayes_R2(norin.brm2, re.form=~(1|FISHID), summary=FALSE) %>%
median_hdci
bayes_R2(norin.brm2, re.form=~(TRIAL|FISHID), summary=FALSE) %>%
median_hdci
emmeans(norin.brm2, pairwise~TRIAL)
norin.em <- norin.brm2 %>%
emmeans(~TRIAL) %>%
pairs() %>%
gather_emmeans_draws() %>%
mutate(Fit=.value)
norin.em %>%
group_by(contrast) %>%
median_hdi()
norin.em %>%
ggplot() +
geom_vline(xintercept=0, linetype='dashed') +
stat_slab(aes(x=.value, y=contrast,
fill = stat(ggdist::cut_cdf_qi(cdf,
.width = c(0.5, 0.8, 0.95),
labels = scales::percent_format())
)), color='black') +
scale_fill_brewer('Interval', direction = -1, na.translate = FALSE) +
theme_bw()
norin.em %>%
group_by(contrast) %>%
summarize(P=sum(Fit>0)/n())
norin.grid <- with(norin, list(SMR_contr=c(min(SMR_contr),
mean(SMR_contr),
max(SMR_contr))))
norin.em <- norin.brm2 %>%
emmeans(~TRIAL|SMR_contr, at=norin.grid) %>%
pairs() %>%
gather_emmeans_draws()
norin.em %>% head
norin.em %>%
group_by(contrast, SMR_contr) %>%
median_hdi()
norin.em %>%
group_by(contrast, SMR_contr) %>%
summarize(P=sum(.value>0)/n())
norin.grid <- with(norin, list(SMR_contr=modelr::seq_range(SMR_contr, n=100)))
newdata <- norin.brm2 %>%
emmeans(~SMR_contr|TRIAL, at=norin.grid) %>%
as.data.frame
head(newdata)
partial.obs <- norin %>%
mutate(Pred = predict(norin.brm2, re.form = NA, summary=TRUE)[,'Estimate'],
Resid = resid(norin.brm2)[,'Estimate'],
Obs = Pred + Resid)
ggplot(newdata, aes(y=emmean, x=SMR_contr, color=TRIAL)) +
geom_point(data=partial.obs, aes(y=Obs)) +
##geom_point(data=partial.obs, aes(y=CHANGE), shape=2) +
geom_ribbon(aes(ymin=lower.HPD, ymax=upper.HPD, fill=TRIAL), alpha=0.3,color=NA) +
geom_line()